%ylevela=xlsread('x:\Eelco Zandberg\Y (in niveaus).xls'); % data in levels
ylevelb=xlsread('x:\Eelco Zandberg\Y (in niveaus minus item 6).xls'); % minus item 6
x=xlsread('x:\Eelco Zandberg\X.xls'); 
WA=xlsread('x:\Eelco Zandberg\Weights matrix (distance capitals).xls'); 
%eig(WA)
% Three W matrices
%WB=xlsread('x:\Eelco Zandberg\Weights matrix (bilateral trade flows).xls');
%WCC=xlsread('x:\Eelco Zandberg\Regional Leaders.xls');
%Wdist=xlsread('x:\Eelco Zandberg\distancematrix.xls');
T=30;
N=62;
%WC=zeros(N,N);
%for i=1:N
%    if (WCC(i,2)~=0 ) WC(i,WCC(i,2))=1; end
%    if (WCC(i,3)~=0 ) WC(i,WCC(i,3))=1; end
%end
%WC=normw(WC);
%Winvdist=1./Wdist;
%for i=1:N
%    Winvdist(i,i)=0;
%    Wdist(i,i)=0;
%end
%Wtwee=normw(Winvdist^2);
%
% Choice of model
%
% Choose one of the two dependent variables 
%ylevel=ylevela;fprintf(1,'Y alle items \n');
ylevel=ylevelb;fprintf(1,'Y minus item 6 \n');
ylevel=100*ylevel;
% Choose one of the W matrices
W=WA;fprintf(1,'W obv distance capitals \n');p=eig(W);eigmax=p(2);peig=1;
%W=Wtwee;fprintf(1,'W obv distance capitals squared\n');p=eig(W);eigmax=p(2);peig=1;
%W=WB;fprintf(1,'W obv bilateral trade flows \n');p=eig(W);eigmax=p(3);peig=1;
%W=WC;fprintf(1,'W obv regionals leaders \n');p=eig(W);eigmax=0;peig=1;
%W=WA5;fprintf(1,'W obv distance capitals < 5,000 km\n');p=eig(W);eigmax=max(p);
%W=WA10;fprintf(1,'W obv distance capitals < 10,000 km\n');p=eig(W);eigmax=max(p);
%W=WA15;fprintf(1,'W obv distance capitals < 15,000 km\n');p=eig(W);eigmax=max(p);

% Choose one of the different sets of control variables
x1=x(:,1:5);
x2=x(:,1:7);
x3=x(:,[1:7,9:12]);
vnames1=strvcat('FL(t)','FL(t-1)','W*FL(t-1)','BANK','CUR','DEBT','RECESSION','HINFL');
vnames2=strvcat('FL(t)','FL(t-1)','W*FL(t-1)','BANK','CUR','DEBT','RECESSION','HINFL','FIRSTYEAR','IMF');
vnames3=strvcat('FL(t)','FL(t-1)','W*FL(t-1)','BANK','CUR','DEBT','RECESSION','HINFL',...
               'FIRSTYEAR','IMF','LEFT','RIGHT','DEMO','GOVFRAC');
%x=x1;vnames=vnames1;fprintf(1,'x var 1-5 \n');
%x=x2;vnames=vnames2;fprintf(1,'x var 1-7 \n');
x=x3;vnames=vnames3;fprintf(1,'x var 1-7 + 9-12 \n');
%
% End choice of model
%
for t=1:T+1
    t1=1+(t-1)*N;t2=t*N;
    Wylevel(t1:t2,1)=W*ylevel(t1:t2,1);
end
y=ylevel(N+1:end)-ylevel(1:end-N);
%
% First the model without time-period fixed effects
%
info.lflag=0;
info.tl=1;
info.stl=1;
info.ted=1; %transformed approach
info.dyn=1;
info.model=1;
results=sar_panel_FE(ylevel(N+1:end),[ylevel(1:end-N) Wylevel(1:end-N) x],W,T,info);
prt_spnew(results,vnames,1);
results1=sar_jihai(ylevel,x,W,info);
results.beta=results1.theta1(1:end-2);
results.rho=results1.theta1(end-1);
results.tstat=results1.tstat(1:end-1);
results.sige=results1.theta1(end); %%%
results.lik = f2_sarpanel(results1.theta1,results1.yt,results1.zt,W,results.lndet,T); %%%
prt_spnew(results,vnames,1);

% Needed for F-test of time-period fixed effects
RRSS=results.sige*N*T;
%
% Now model with time-period fixed effect
%
info.lflag=0;
info.tl=1;
info.stl=1;
info.ted=1; %transformed approach
info.dyn=1;
info.model=3;
results=sar_panel_FE(ylevel(N+1:end),[ylevel(1:end-N) Wylevel(1:end-N) x],W,T,info);
%prt_spnew(results,vnames,1);
results1=sar_jihai_time(ylevel,x,W,info);
results.beta=results1.theta1(1:end-2);
results.rho=results1.theta1(end-1);
results.tstat=results1.tstat(1:end-1);
results.sige=results1.theta1(end); %%%
results.lik = f2_sarpanel(results1.theta1,results1.yt,results1.zt,W,results.lndet,T-1); %%%
prt_spnew(results,vnames,1);
varcov=results1.varcov;
R=zeros(1,1);
npar=length(results1.theta1);
tau = results1.theta1(1,1);
eta = results1.theta1(2,1);
rho = results1.theta1(npar-1,1);
R(1)=tau+rho+eta-1;
Rafg=zeros(1,npar);
Rafg(1,1)=1;Rafg(1,2)=1;Rafg(1,npar-1)=1;
sumpar=tau+rho+eta;
Wald=R(1)'*inv(Rafg(1,:)*varcov*Rafg(1,:)')*R(1);
F1=1-chis_cdf (Wald,1);
[sumpar Wald F1]
% F-test of time-period fixed effects
URSS=results.sige*N*T;
[Kjunck K1]=size([ylevel(1:end-N) Wylevel(1:end-N) x Wylevel(N+1:end)]);
F0=((RRSS-URSS)/(T-1))/(URSS/((N-1)*(T-1)-K1))
kansfo = fdis_prb(F0, T-1, (N-1)*(T-1)-K1)
%
% Spatial first-differences
%
sigma=(eye(N)-W)*(eye(N)-W)';
[V D]=eig(sigma);
transf=D(1+peig:end,1+peig:end)^(-0.5)*V(:,1+peig:end)'*(eye(N)-W);
Wstar=D(1+peig:end,1+peig:end)^(-0.5)*V(:,1+peig:end)'*W*V(:,1+peig:end)*D(1+peig:end,1+peig:end)^(0.5);
info.lflag=0;
for t=1:T+1
    t1=1+(t-1)*N;t2=t*N;
    ts1=1+(t-1)*(N-peig);ts2=t*(N-peig);
    ystar(ts1:ts2,1)=transf*ylevel(t1:t2,1);
    Wystar(ts1:ts2,1)=Wstar*ystar(ts1:ts2,1);
end

for t=1:T
    t1=1+(t-1)*N;t2=t*N;
    ts1=1+(t-1)*(N-peig);ts2=t*(N-peig);
    xstar(ts1:ts2,:)=transf*x(t1:t2,:);
end
info.lflag=0;
info.tl=1;
info.stl=1;
info.dyn=1;
info.model=1;
N1=N-peig;
results=sar_panel_FE(ystar(N1+1:end),[ystar(1:end-N1) Wystar(1:end-N1) xstar],Wstar,T,info);
%prt_spnew(results,vnames,1);
results1=sar_jihai(ystar,xstar,Wstar,info);
results.beta=results1.theta1(1:end-2);
results.rho=results1.theta1(end-1);
results.tstat=results1.tstat(1:end-1);
results.sige=results1.theta1(end); %%%
results.lik = f2_sarpanel(results1.theta1,results1.yt,results1.zt,Wstar,results.lndet,T-1); %%%
prt_spnew(results,vnames,1);

tau = results1.theta1(1,1);
eta = results1.theta1(2,1);
rho = results1.theta1(npar-1,1);
varcov=results1.varcov;
btemp=results1.theta1;
R=zeros(1,1);
R(1)=tau+rho+eta-1;
Rafg=zeros(1,npar);
Rafg(1,1)=1;Rafg(1,2)=1;Rafg(1,npar-1)=1;
sumpar=tau+rho+eta;
check=tau+eigmax*(rho+eta)-1; % this should be negative
Wald=R(1)'*inv(Rafg(1,:)*varcov*Rafg(1,:)')*R(1);
F1=1-chis_cdf (Wald,1);
[sumpar check Wald F1]
%
% Calculation of direct/indirect effects
%
NSIM=1000;
simresults=zeros(npar-2,NSIM);
simdir=zeros(npar-3,NSIM);
simind=zeros(npar-3,NSIM);
simtot=zeros(npar-3,NSIM);
for sim=1:NSIM
    parms = chol(real(varcov)+0.001)'*randn(size(btemp)) + btemp;
    rhosim = parms(npar-1,1);
    betasim = parms(3:npar-2,1);
    etasim=parms(1,1);
    tausim = parms(2,1);
    simresults(:,sim)=[tausim;betasim;rhosim];
    SC=inv(eye(N)-rhosim*W)*((tausim-1)*eye(N)+(rhosim+etasim)*W);
    p=1;
    EAVD(p,1)=sum(diag(SC))/N; % average direct effect
    EAVI(p,1)=sum(sum(SC,2)-diag(SC))/N; % average indirect effect
    EAVC(p,1)=sum(sum(SC,1)'-diag(SC))/N; % average indirect effect
    simdir(p,sim)=EAVD(p,1);
    simind(p,sim)=EAVI(p,1);
    simtot(p,sim)=EAVD(p,1)+EAVI(p,1);
    for p=2:npar-3
        C=zeros(N,N);
        for i=1:N
            for j=1:N
                if (i==j) C(i,j)=betasim(p-1);
                end
            end
        end
        S=inv(eye(N)-rhosim*W)*C;
        EAVD(p,1)=sum(diag(S))/N; % average direct effect
        EAVI(p,1)=sum(sum(S,2)-diag(S))/N; % average indirect effect
        EAVC(p,1)=sum(sum(S,1)'-diag(S))/N; % average indirect effect
        simdir(p,sim)=EAVD(p,1);
        simind(p,sim)=EAVI(p,1);
        simtot(p,sim)=EAVD(p,1)+EAVI(p,1);
    end
end
[mean(simresults,2) mean(simresults,2)./std(simresults,0,2)]
fprintf(1,'First effect is convergence effect \n');
[mean(simdir,2) mean(simdir,2)./std(simdir,0,2) mean(simind,2) mean(simind,2)./std(simind,0,2)...
    mean(simtot,2) mean(simtot,2)./std(simtot,0,2)]
%
% Impose restriction
%
theta2=btemp-varcov*Rafg'*inv(Rafg(1,:)*varcov*Rafg(1,:)')*R(1);
varcov2=varcov-varcov*Rafg'*inv(Rafg(1,:)*varcov*Rafg(1,:)')*Rafg*varcov;
results2=results;
results2.meth='sar_jihai_restricted';
results2.theta1=theta2;
results2.tstat1=theta2./(sqrt(abs(diag(varcov2))));
results2.sige=theta2(end); %%%
results2.lik = f2_sarpanel(results2.theta1,results1.yt,results1.zt,Wstar,results.lndet,T); %%%
help=theta2(end-1)*kron(speye(T),Wstar)*results1.yt;
residr=results1.yt-help-results1.zt*theta2(1:end-2);
yme=results1.yt-mean(results1.yt);
rsqr2=yme'*yme;
rsqr1 = residr'*residr;
results2.rsqr=1.0-rsqr1/rsqr2; %rsquared
res1=yme;
res2=((speye((N1)*T))-theta2(end-1)*kron(speye(T),Wstar))\(results1.zt*theta2(1:end-2))-mean(results1.yt);
rsq1=res1'*res2;
rsq2=res1'*res1;
rsq3=res2'*res2;
results2.corr2=rsq1^2/(rsq2*rsq3); %corr2
prt_sardynamic(results2,vnames,1);
tau = results2.theta1(1,1);
eta = results2.theta1(2,1);
rho = results2.theta1(npar-1,1);
%
% Calculation of direct/indirect effects
%
clear SC;
NSIM=1000;
simresults=zeros(npar-2,NSIM);
simdir=zeros(npar-3,NSIM);
simind=zeros(npar-3,NSIM);
simtot=zeros(npar-3,NSIM);
for sim=1:NSIM
    parms = chol(real(varcov2)+0.001)'*randn(size(theta2)) + theta2;
    rhosim = parms(npar-1,1);
    betasim = parms(3:npar-2,1);
    tausim = parms(2,1);
    simresults(:,sim)=[tausim;betasim;rhosim];
    SC=(tausim-1)*inv(eye(N)-rhosim*W)*(eye(N)-W);
    p=1;
    EAVD(p,1)=sum(diag(SC))/N; % average direct effect
    EAVI(p,1)=sum(sum(SC,2)-diag(SC))/N; % average indirect effect
    EAVC(p,1)=sum(sum(SC,1)'-diag(SC))/N; % average indirect effect
    simdir(p,sim)=EAVD(p,1);
    simind(p,sim)=EAVI(p,1);
    simtot(p,sim)=EAVD(p,1)+EAVI(p,1);
    for p=2:npar-3
        C=zeros(N,N);
        for i=1:N
            for j=1:N
                if (i==j) C(i,j)=betasim(p-1);
                end
            end
        end
        S=inv(eye(N)-rhosim*W)*C;
        EAVD(p,1)=sum(diag(S))/N; % average direct effect
        EAVI(p,1)=sum(sum(S,2)-diag(S))/N; % average indirect effect
        EAVC(p,1)=sum(sum(S,1)'-diag(S))/N; % average indirect effect
        simdir(p,sim)=EAVD(p,1);
        simind(p,sim)=EAVI(p,1);
        simtot(p,sim)=EAVD(p,1)+EAVI(p,1);
    end
end
[mean(simresults,2) mean(simresults,2)./std(simresults,0,2)]
fprintf(1,'First effect is convergence effect \n');
[mean(simdir,2) mean(simdir,2)./std(simdir,0,2) mean(simind,2) mean(simind,2)./std(simind,0,2)...
    mean(simtot,2) mean(simtot,2)./std(simtot,0,2)]
clear simresults simdir simind simtot;
%end